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1. Introduction 

Neutrinos are among the most abundant particles in our Universe and therefore play an 
important role in the formation of large-scale structure. In the past 10 years neutrino 
oscillation experiments have verified that neutrinos have small, but non-zero masses and 
that they therefore contribute to the dark matter density in the Universe. 

Neutrino oscillation experiments have established the existence of two distinct 
neutrino mass splittings, Am 2 ,- = mf — raj. Am 2 2 ~ 7 x 10~ 5 eV 2 has been found 
by observing Solar and reactor neutrinos, and |Am| 3 | ~ 2.5 x 10 _3 eV 2 by observing 
atmospheric neutrinos and by accelerator experiments (see e.g. [U [2] for recent data 
analyses). This gives a lower bound on the sum of the neutrino masses which is 
^2 m v > 0.05 eV for the normal hierarchy and ^2 m v > 0.1 eV for the inverted hierarchy. 

Using a combination of cosmological observables such as the cosmic microwave 
background, large scale structure, and type la supernovae a bound on the sum of 
neutrino masses of ^ 0-2 — 0.7 eV (95% C.L.) has been derived [31 IH El El 13 

El El CD]- The actual value of the bound depends both on the specific combination of 
data sets used and on the cosmological model. At the moment a fairly robust upper 
limit can be taken to be roughly m v < 0.6 — 0.7 eV. 

This bound mainly uses structure formation data in the linear regime (k < 
0.15 hMpc~ l ), but even on such large scales some non-linear contamination is present 
and must be modelled. This has been done in detail for variations of the standard 
ACDM model, but so far not in a truly consistent way for models which include massive 
neutrinos. 

For the study of neutrino mass bounds with present data the current precision of 
theoretical power spectrum calculations in the regime k ~ 0.1 — 1 /iMpc -1 is sufficient 
in the sense that the final result does not depend very significantly on the method 
used (however, see e.g. pj] for an exception). However, future high precision lensing 
and galaxy redshift surveys, such as the LSST [12], will constrain the matter power 
spectrum to percent level precision on scales of k ~ 0.1 — 1 /iMpc -1 . 

At this level of accuracy the inclusion of neutrinos in the calculation of the matter 
power spectrum will be crucial, even for neutrino masses close to the current lower 
bound. While the effect of neutrinos in the linear theory power spectrum is understood 
in detail much work remains to be done in the semi-linear and non-linear regimes. 
Furthermore, semi-analytic models of non-linear structure formation such as the halo 
model must always be checked against iV-body simulations. 

Neutrinos have thermal velocities which exceed typical gravitational flow velocities 
by an order of magnitude or more, even at low redshift. Therefore, when setting up 
initial conditions for A-body simulations the neutrino thermal velocities need to be 
added to the gravitational flow velocities found from the neutrino transfer function. 
Including neutrino thermal velocities in A-body simulations is not without problems 
because the thermal motion introduces a noise term which will affect the underlying 
gravitational evolution unless high resolution is used. 



The Effect of Thermal Neutrino Motion on the Matter Power Spectrum 



3 



The problem is similar to the well known problem that particles in an iV-body 
simulation must be put in a configuration like a grid or a glass which respects small- 
scale homogeneity in order not to introduce a large Poisson noise term. When 
thermal velocities are included momentum phase space must be sampled and the initial 
thermal velocities should preferably be chosen in a way where small-scale momentum 
conservation is enforced. 

In the present paper we present the first A^-body simulations in which neutrino 
thermal velocities have been included and convergent results achieved for k < 
1.5/iMpc -1 . In the early 90s [TBI H3] included neutrino thermal velocities albeit for 
a model with much higher neutrino masses, and no convergent results were proven. We 
also note another recent paper [TH] in which thermal motion has been included for warm 
dark matter, a set-up somewhat different from what is studied here. Neutrino clustering 
in the small-scale limit has been studied using a completely different approach [161 E] 
where the full Boltzmann equation has been solved for neutrinos in a background given 
by pure cold dark matter (CDM), i.e. no feed-back is included and neutrinos are treated 
as tracer particles in the underlying gravitational potential. This method allows for 
arbitrarily high resolution, but is not truly self-consistent. 

In Section [2] we describe the linear theory evolution and the set-up of initial 
conditions for our neutrino iV-body simulations. In Section [3] we describe how our 
iV-body simulations are performed, and in Section H] the results are described in detail. 
Finally, Section [5] contains a discussion and our conclusions. 

2. Linear Evolution of Perturbations and Initial Conditions 

2.1. Linear Theory 

Evolving the primordial density perturbations set down by inflation to the present 
involves several steps. As long as perturbations remain small the evolution can be 
calculated precisely using the linearised Einstein and Boltzmann equations [18], using 
software such as CAMB [22] or CMBFAST [23]- 

In this regime the power spectrum can be factorised into a primordial component, 
Po{k), and a transfer function (TF), T(k,z), which contains all information about the 
evolution of structure so that P(k, z) = Po(k)T 2 (k, z). 

However, once structure enters the non-linear regime precision studies require the 
use of A^-body simulations. To set up the initial conditions (ICs) for our simulations we 
have calculated the TFs using CMBFAST (using CAMB yields similar results). 

We have assumed a flat cosmological model with density parameters fi& = 0.05, 
Q m = 0.30 and Q\ = 0.70 for the baryon, matter and cosmological constant components, 
respectively, and a Hubble parameter of h = 0.70. We vary the CDM and neutrino 
density parameters (ficDM and fi u , respectively) fulfilling the condition I^cdm + Q u = 
0.25. We have assumed a primordial power spectrum of the form Po(k) = Ak n with 
n — 1, i.e. a standard scale-invariant Harrison-Zel'dovich spectrum. The amplitude, 
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Figure 1. The linear theory transfer functions at z = 4 for the CDM and neutrino 
components. 



giving erg = 0.89 for a pure ACDM model, is chosen so as to fit the CMBR [26] on large 
scales. 

Different TFs at a redshift of 4 are shown in Fig. [TJ The effect of neutrino free- 
streaming on the TF is clearly seen, and is more pronounced in the lower mass neutrino 
case. 

2.2. Initial Conditions with Two Species 

The TFs are used to generate the position and velocity ICs for the A-body particles. 
The Zel'dovich Approximation (ZA) [29], based on first-order Lagrangian perturbation 
theory, is a standard way to calculate the ICs. Using second-order Lagrangian 
perturbation theory (2LPT) [191 ED] it is possible to generate ICs in the quasi-linear 
regime. We are particularly interested in these second-order corrections because the 
presence of neutrino thermal velocity noise requires simulations to be started at low 
redshift, close to the non-linear regime. 

We have used a modified version of the 2LPT initial conditions code of [21] to 
generate the ICs for our simulations. The modifications are described in detail below. 

In the A"-body simulations we include two particle species, CDM and neutrinos. 
It is not possible to generate particle ICs for each species simultaneously with the 
ZA+2LPT formalism. Instead we generate ICs one particle species at a time by using 
their respective TFs (for CDM we have used a weighted sum of the CDM and baryon 
TFs for consistency, even though this only matters at high redshift). Thereafter, the 
A-body particle masses for each species are scaled so that they make up their proper 
fraction of Q m . 

The initial positions for the A-body particles are found by adding a displacement 
field to a regular grid. A standard way to get the particle velocities is to differentiate 
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Figure 2. Images of the CDM and neutrino density distributions in a slice of the 
simulation volume. The images span 512 h~ 1 Mpc on a side and has a depth of 
10ft. _1 Mpc. To produce the images we have interpolated the masses of the TV-body 
particles to a regular grid with the adaptive smoothing kernel of [25] . The images show 
the densities for the CDM component (left), neutrinos with X) m ^ = 0.6 eV (middle), 
and neutrinos with ^ m„ = 0.3 eV (right). The top row is at Zi = 4 and the bottom 
row at z = 0. To enhance the dynamic range of the CDM structures the square root 
has been taken of the CDM density field in the z = image. The Yl m v = 0-3 eV 
neutrino image at z = displays artificial small-scale structures in the voids caused by 
neutrino A-body particle shot-noise. All the images are made from simulations with 
512 3 neutrino A-body particles. 

the initial position displacements. This procedure involves using several numerically 
determined fitting factors, and therefore breaks down when two species with different 
TFs are present since then the growth factor is both species and mode dependent. 
Instead, we get the velocities by generating two displacement fields centered around our 
starting redshift and then take the time difference. We have tested that these velocities 
do not depend on the distance in redshift between the two extra displacement fields in 
a suitable range around our starting redshift. 

2LPT involves a relation between the first- and second-order growth factors. But 
since the perturbed energy density even at = 4 (zi designates the A^-body starting 
redshift) is vastly dominated by CDM, we can neglect the neutrino contribution to the 
driving term for the CDM growth factor since this would give a small correction to a 
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second-order term. Because the neutrinos are smoothly distributed, first-order theory 
is very accurate for them, and therefore the second-order corrections give negligible 
contributions to the neutrinos. We note that this is somewhat similar to the approach 
taken in Refs. [101 (23> [2H] for the case of semi-analytic power spectrum modelling. 

In the top row of Fig. [2] the CDM and neutrino distributions are shown at Zi = 4. 
Since we run simulations with different numbers of CDM and neutrino iV-body particles, 
the amplitude and phase of the common large-scale modes of each species have been 
generated with the same random numbers. Therefore, from Fig. [2] one can identify the 
components as having the same large-scale structure. Because of the higher value of 
the CDM TF at small scales, this component has additional clearly visible small-scale 
structure. As expected, the neutrino distributions are more homogeneous than their 
CDM counterpart. The Yl m v = 0.6 eV and = 0.3 eV distributions are similar at 

Zi = 4, with the latter distribution being marginally more homogeneous, as can also be 
inferred from the TFs in Fig. [TJ 



2.3. Thermal Velocities 

Since the neutrino mass is small the thermal velocities cannot be neglected even at low 
redshift if percent level precision is to be obtained. Instead, the velocity of a given 
neutrino iV-body particle contains a thermal component drawn from an equilibrium 
Fermi-Dirac (FD) distribution. Assuming isotropy for the thermal component the 
probability, Pr, for a neutrino having a momentum smaller than p is given by 

Pr(< p)=N -— ^ -dp', (1) 

where N is a normalisation which ensures that the probability is bounded between 
and 1. Fig. [3] shows the cumulative FD distributions for different one-particle neutrino 
masses. The magnitudes of the neutrino thermal velocities were drawn randomly from 
the FD velocity distribution, and the directions of the thermal velocities were drawn at 
random as well. This reflects the fact that the neutrino velocity can be split into two 
uncorrelated contributions. An equilibrium, random contribution arising from the FD 
distribution and an out-of-equilibrium gravitational flow contribution found from the 
ZA+2LPT formalism. 

The neutrinos decoupled from the baryon-photon plasma before the annihilation 
of electrons and positrons, and as a result the neutrino temperature is related to the 
photon temperature by the approximate relation T u ~ T 7 (4/ll) 1//3 . Assuming p <C m v 
at late times the 3 neutrino generations contribute the following to the present neutrino 
density parameter [301 EH E21 [33l EH ESI EEl [37] 

= 93.8tf eV (2) 
The finite number of neutrino TV-body particles gives rise to noise in the power 
spectrum. Physically this leads to spuriously strong clustering of neutrinos on small 



The Effect of Thermal Neutrino Motion on the Matter Power Spectrum 7 




v [10 3 km s" 1 ] 

Figure 3. The cumulative Fermi-Dirac distributions as a function of velocity for most 
of the simulations listed in Table [1] 

scales, even at early times. As will be discussed below we have carefully checked that 
this noise term is negligible on the scales of interest in our simulations. 

However, since the neutrino thermal velocity distribution redshifts as a" 1 , and the 
neutrino and CDM out-of-equilibrium density contrasts and velocity fields increase with 
the expansion of the Universe, it is paramount to start the iV-body simulations as late, 
and therefore as close to the non-linear regime, as possible. This in turn constrains the 
scales on which clustering can be reliably calculated. 

We have tested the validity of starting the TV-body simulations as late as Zi = 4. 
This has been done by starting 4 pure ACDM simulations, at Zi = 49, 11.5, 7 and 4, and 
evolving them to the present. The difference in the matter power spectrum between these 
simulations is at the few percent level on small scales. Since the neutrino simulations 
include a smaller amount of CDM this is an upper limit. The 2LPT corrections were 
crucial for achieving such a small discrepancy. Since we are interested in quantifying 
the relative effect of including neutrinos, and not the exact absolute value of the matter 
power spectrum, we are justified in choosing z^ = 4 as the lowest starting redshift for 
our TV-body simulations. If the absolute power spectrum is desired at the percent level 
a higher starting redshift should be chosen, but we stress that this will have a minimal 
effect on the relative change in the matter power spectrum coming from including the 
neutrino component. 

3. iV-body Simulations 

The TV-body simulations were performed with the publicly available TV-body code 
gadget- 2 [2T] run in the hybrid TreePM mode. Gas physics was neglected, since 
it does not significantly affect the scales of interest. 

The TV-body simulations are evolved with Newtonian dynamics. There are several 
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qO cO cO cO 
°1 J 2 J 3 J 4 


cO.15 cO.15 cO.15 cO.15 cO.15 cO.15 
J l J 2 °3 °i °5 °6 


c0.3 cO.3 cO.3 cO.3 cO.3 
°1 D 2 D 3 °4 J 5 


qO.45 cO.45 cO.45 
°1 °2 °3 




Y.'mv [eV] 

n„,o [%] 

FD 





No No No No 
4 7 11.5 49 


256 3 128 3 256 3 256 3 
0.15 0.15 0.15 0.15 0.15 0.15 
0.33 0.33 0.33 0.33 0.33 0.33 

No No Yes Yes No No 
4 4 4 4 49 49 


256 3 128 3 256 3 512 3 
0.3 0.3 0.3 0.3 0.3 
0.65 0.65 0.65 0.65 0.65 

No No Yes Yes Yes 
4 4 4 4 4 


256 3 256 3 
0.45 0.45 0.45 
0.98 0.98 0.98 

No No Yes 
4 4 4 




c0.6 cO.6 c/0.6 cO.6 cO.6 cO.6 cO.6 cO.6 
°1 °2 °3 °4 °5 °6 °7 °S 


qO.6 qO.6 cO.6 cO.6 
°9 °10 J ll J 12 


c-0.6 c-0.6 cO.6 cO.6 
°13 °U °15 J 16 


cO.6 c-0.6 c/0.6 
°17 °18 °19 


[eV] 

n„,o [%] 

FD 


256 3 32 3 64 3 128 3 256 3 512 3 1024 3 
0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 
1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 

No No Yes Yes Yes Yes Yes Yes 
4444444 4 


256 3 128 3 256 3 
0.6 0.6 0.6 0.6 
1.3 1.3 1.3 1.3 
No No Yes Yes 
7 7 7 7 


256 3 128 3 256 3 
0.6 0.6 0.6 0.6 

I. 3 1.3 1.3 1.3 
No No Yes Yes 

II. 5 11.5 11.5 11.5 


128 3 256 3 512 3 
0.6 0.6 0.6 
1.3 1.3 1.3 

Yes Yes Yes 
49 49 49 



Table 1. Parameters for our iV-body simulations. N v is the number of neutrino iV- 
body particles. ^ m v is the total neutrino mass, and it is in all cases related to the 
one-particle neutrino mass, m„, by — im u . fl U fl is the fraction of the critical 

density contributed by the neutrinos today (see Eq. [2]), and FD indicates whether or 
not the neutrinos have been given a thermal velocity from the relativistic Fermi-Dirac 
distribution (see Eq. Q]). Finally, %{ indicates the A^-body starting redshift. Note that 
N u = combined with ^2 m v > indicates that the neutrino distribution has been 
kept totally homogeneous in the iV-body simulation. 

reasons for this. First, the ZA+2LPT formalism is stricly Newtonian. Second, even 
though the neutrino thermal velocities approach up to 15% of the speed of light, the 
relativistic corrections are much smaller than the desired accuracy. 

The simulations are listed in Table [H and include 256 3 CDM particles in a 
512/i -1 Mpc box. By running a simulation with 512 3 CDM particles we have tested 
the validity of using 256 3 CDM particles. The simulation with the highest number of 
iV-body particles gave a correction to the matter power spectrum which was at the one 
percent level on small scales. Since we only quantify the relative effect of including 
neutrinos, 256 3 CDM particles is sufficient. To test the convergence of the matter and 
neutrino power spectra we have made runs with different numbers of neutrino iV-body 
particles. 

As a reference, and in order to test how the neutrino thermal velocity affect the 
power spectra, we have run simulations with 256 3 neutrinos with no thermal velocities 
included. 

We have also tested how much the neutrinos contribute to the matter power 
spectrum by running simulations with the neutrinos included in the linear evolution 
but neglecting the perturbed neutrino component in the iV-body simulations. Here, the 
neutrinos were still included in the calculation of the Hubble parameter so that the 
background evolution is identical. Physically, this corresponds to keeping the neutrino 
distribution totally homogeneous in the iV-body simulation. 

The gravitational Tree part in gadget-2 sets the timestep according to the 
gravitational acceleration, not the neutrino thermal velocity, so that the fast moving 
neutrinos may not be accurately evolved. We have tested if the neutrino timestep was 
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small enough to simulate the scales of interest by decreasing it by roughly a factor of 5 
in the Tree part. The correction found was well below 0.1%. This very small correction 
can be explained by the fact that the neutrino small-scale structure does not contribute 
to the matter power spectrum on the smallest scales simulated. 

Finally, we have investigated the effect of the finite box size for the case with 
^2m u = 0.6 eV and Zj = 4. This has been done by reducing the size of the simulation 
volume as well as the number of CDM particles by a factor of two. A different set of 
random numbers was used to generate the initial conditions. The amount of suppresion 
due to neutrinos and the position of the turnover in the difference power spectrum were 
found to be consistent with the results presented in Fig. HI 

4. Results 

4-1. Damping and Convergence of the Power Spectrum 

Fig. H] shows the well-known scale dependent damping of power when neutrinos are 
included. As expected, the damping is much larger in the ~Y^m v = 0.6 eV case as 
compared to ^2m u = 0.15 eV because ^cdm is smaller in the former case even though 
the neutrino thermal velocity is largest for the lower mass neutrino. 

From Fig. H]it can be seen that in the Yl m f = 0.6 eV case the damping is almost 
independent of Zj. Note that the relative decrease in power for a given simulation with 
neutrinos is taken with respect to a ACDM simulation started at the same Zj. This 
has been done to remove the dependence on the starting redshift when comparing the 
models. 

Fig. H] also shows the damping expected from linear theory, which is accurate out 
to k ~ 0.2h Mpc" 1 . On smaller scales non-linear theory predicts a substantially larger 
damping of the power spectrum. The departure from linear theory increases with higher 
^2 m u . Note that non-linear theory predicts a turnover in the difference power spectrum. 
At the starting redshift, Zj, this turnover is not present, neither in the TFs nor when 
the ICs have been calculated. 

From a mode-coupling point of view the turnover occurs because the small-scale 
modes (k ~ 1/iMpc -1 ) in the simulations with neutrinos get relatively more out of 
their coupling to the large-scale modes than do the same modes in the pure ACDM 
simulations. The turnover does not appear because of mode-coupling between the 
small-scale modes themselves. Physically, large-scale potential gradients are needed to 
make the small-scale structures grow substantially, and the large-scale perturbations 
are less affected by neutrino free-streaming than are the small-scale perturbations. 
The non-linear collapse of the halos then gives mode-coupling and, combined with 
the characteristic damping of the difference power spectrum caused by neutrino free- 
streaming, gives the turnover. 

The turnover is not caused by clustering of the neutrino component on small scales. 
This can be seen by comparing the two figures in the top panel of Fig. [5] (see below for 
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Figure 4. Relative differences in the matter power spectra at z = between pure 
ACDM models and models with neutrinos included. The differences expected from 
linear theory are also shown. The horizontal black lines indicate a relative power 
spectrum suppression of —9.8 Q, v /£l m . 



further comments on Fig. If small-scale neutrino clustering did contribute to the 
matter power spectrum on scales k ~ 1 h Mpc -1 then the matter power spectrum should 
be affected by the different noise levels seen in the neutrino power spectra. But since 
the matter power spectrum has converged for 128 3 to 1024 3 neutrino iV-body particles 
the turnover is not caused by neutrino small-scale clustering. 

The characteristic scale at which the turnover appears, fc turn , moves to smaller k as 
Q u is decreased. k tnTn is not directly proportional to fl u and therefore to the neutrino 
free-streaming length. Instead, k tUTn is related to the amount of ^cdm- Since we keep 
Q m fixed, decreasing Q u increases the amount of CDM, and more CDM makes the 
halos collapse on larger scales. Since the turnover is related to the non-linear collapse 
of structures, /c tU rn moves to larger scales as VL V is decreased. Since there is a linear 
boundary condition at k ~ 0.2h Mpc -1 , the turnover cannot propagate beyond this 
value for fl u approaching zero. 

The maximum relative magnitude of the power spectrum suppression is roughly 
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Figure 5. Top left: Differences in % in the matter power spectra at z = with 
^2 rn v = 0.6 eV neutrinos and z\ = 4. The differences are taken with respect to 
the 1024 3 neutrino simulation. Top right: Neutrino power spectra at z = with 
m„ = 0.6 eV neutrinos and Zj = 4. Bottom left: Differences in % in the matter 
power spectra at z — with ^ m„ = 0.3 eV neutrinos and z% = 4. The differences 
are taken with respect to the 512 3 neutrino simulation. Bottom right: Neutrino power 
spectra at z = with ^ m„ = 0.3 eV neutrinos and z% = 4. 



given by 

AP n u 

~p ~ " 9 - 8 ?T' (3) 

max J 

which is about 20% larger than the linear theory prediction of AP/P\\ [n ~ — 8Q U /Q m 
With future high precision measurements of the matter power spectrum 
approaching the 1% precision on these scales this effect must be taken into account 
even for hierarchical neutrino masses. 

The differences for the matter power spectra in the m v = 0.6 eV function 
of different numbers of neutrino iV-body particles are shown in the top left panel of 
Fig. El By using different mass assignment methods, Nearest-Grid-Point, deconvolved 
Cloud- In- Cell, the adaptive smoothing kernel of [25J as well as dividing the simulation 
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Figure 6. The effect on the matter power spectrum at z = of neglecting thermal 
velocities in the iV-body simulations. The neutrinos have been correctly included in 
the linear evolution. 



volume into equal cubical cells and superposing the particle distributions for attaining 
higher spacial resolution, we have tested that the power spectra do not have any artificial 
features caused by the particular mass assignment method used. For all neutrino masses 
simulated it is necessary to include at least 128 3 neutrinos in the iV-body simulation and 
give them a thermal velocity to calculate the matter power spectrum at percent level 
accuracy on the relevant scales. For the m u = 0.6 eV case and Zj = 49 it is furthermore 
necessary to include 512 3 neutrino iV-body particles to achieve convergence. 

The neutrino power spectra for ^ m v = 0.6 eV as a function of different numbers 
of neutrino iV-body particles are shown in the top right panel of Fig. [5l On small 
scales the neutrino power spectrum is in most cases completely flat because the coarse 
sampling of the neutrino velocity distribution introduces a white noise term. That this 
is indeed the reason can be seen from the fact that the noise level decreases rapidly 
when more neutrino iV-body particles are used. For runs with 128 3 to 1024 3 neutrino 
vY-body particles the noise effect on the matter power spectrum is kept well below the 
1% level on relevant scales. It can also be seen that the neutrino power spectra for 
N v = 256 3 and N v = 512 3 converge out to k ~ 0.2 h Mpc" 1 and that the N v = 512 3 and 
N u = 1024 3 simulations converge out to k ~ 0.5 h Mpc" . 

In the bottom panel of Fig. [3] the corresponding convergence of the matter power 
spectra (left) and the neutrino power spectra (right) for the ^2 m v = 0-3 eV case are 
shown. The same general trends from the m v = 0.6 eV case can be seen. 

4-2. The Effect of Neglecting the Thermal Component 

Fig. [6] shows the effect of neglecting the neutrino thermal velocity in the simulations. 
In the figure the power spectra are normalized with respect to the simulation with the 
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Figure 7. Left: The effect of neglecting the perturbed neutrino component at z = 
in the iV-body simulations. Right: The effect of neglecting the thermal velocities or 
the perturbed neutrino component for J2 m v = 0.15 eV and a starting redshift of 49. 
The 'Base' has been estimated, as explained in the text. In all cases the neutrinos have 
been correctly included in the linear evolution. 



highest number of neutrino iV-body particles for a given mass and starting redshift. 
As expected, neglecting the thermal velocity component increases the amplitude of the 
power spectrum, because without thermal velocities the neutrinos act as an extra CDM 
species, and the effect of free-streaming is neglected once the simulation is started. The 
effect is more pronounced at smaller scales, and is largest in the highest mass neutrino 
case, since it contributes a larger fraction of Q m , even though it has a smaller omitted 
thermal velocity component than in the cases of the lower mass neutrinos. 

For the case of tti v = 0.6 eV the effect is as large as ~ 5% for Zi — 4 increasing to 
~ 10% for Zi = 11.5 on scales relevant for future large-scale surveys. For higher starting 
redshifts the effect is even larger. We note that for a given starting redshift the effect is 
proportional to m u as expected. 

4-3. High Z\ Low m v Approximation 

Including neutrinos with a total mass in the range 0.3 ~ 0.6 eV in TV-body simulations 
can be done consistently even at redshifts required for getting the absolute power 
spectrum accurately. However, for very light neutrinos, m v ~ 0.15 eV, the thermal 
velocities are semi-relativistic at z — 49. This will not only render Newtonian dynamics 
inaccurately but will also, depending on the number of neutrino iV-body particles used, 
erase the initial neutrino large-scale structure since typical neutrino gravitational flow 
velocities are of order ~ 5 kms -1 , much smaller than the thermal velocities of individual 
particles. Furthermore, the high neutrino velocities make the ./V-body simulation 
timestep very short, increasing simulation times substantially. 

For these reasons it is desirable to develop an approximate method which can be 
used to calculate the matter power spectrum at the 1% level without including thermal 
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velocities. In Fig. [7J (left) we show the effect of neglecting the perturbed neutrino 
component in the TV-body simulations. The CDM iV-body particle mass is still scaled 
according to the total matter density. As expected this gives less power because of the 
lack of feed-back from the large-scale modes of the neutrino component to the CDM 
component. By increasing Zi the neutrino perturbations are omitted over a longer time 
span and therefore affect and decrease the power spectrum more. 

It is important to notice, that this lack of power is very small in the ^ m v = 0.15 eV 
case for k > 0.1 /iMpc -1 , i.e. on scales where non-linear corrections become important. 

Ideally, the ratio of a pure ACDM power spectrum to that of a power spectrum 
with a given neutrino mass, Pacdm/Pm.cdm, should be independent of iV-body starting 
redshift, as long as z,j is not too low. The presence of non-linearities at the time when the 
ICs are calculated make this power spectrum ratio marginally dependent on z%, which 
is also visible in Fig. H]in the Yl 771 ^ = 0.6 eV case. 

Now, focusing on the J2m u = 0.15 eV case and Zi = 49, we would like to estimate 
the error made by not including the perturbed neutrino component in the iV-body 
simulation. Because the damping of the power spectrum is nearly independent of 
Zi it is possible to estimate the matter power spectrum with neutrinos and thermal 
velocities without actually performing the simulation. More specifically, assuming that 
the damping effect is independent of Zj, the matter power spectrum at z = for an 
iV-body simulation started at Zj = 49 with neutrinos and thermal velocities can be 
found from PacdmO; = 4)/P,a CDM (z; = 4) = PacdmO; = 49)/P, AC dmO; = 49). This 
estimated power spectrum is used as the 'Base' in the right panel of Fig. [7J This 
figure shows the error made by neglecting the perturbed neutrino component. On scales 
relevant for iV-body simulations, i.e. k > 0.1 h Mpc" 1 , the error made is at the 1% level. 
For comparison the effect of neglecting the thermal velocities is also shown, and it can 
be seen that it is a better approximation to neglect the perturbed neutrino component 
than to include it but neglect the thermal velocity component. 

In the Y l m v = 0.6 eV case we have all the relevant power spectra at starting 
redshifts of 4 and 49. Therefore we can test the validity of estimating the power spectrum 
with neutrinos and thermal velocities. We have found that the 'Base' in the right panel 
of Fig. [7J is overestimated so that the error made by neglecting the perturbed neutrino 
component is an upper bound. 

We have tried to decrease the error made by neglecting the perturbed neutrino 
component. This has been done by including the perturbed neutrino component 
without thermal velocities in the iV-body simulation, but switching off the effect of 
the gravitational field on the neutrinos, so that the gravitational flow velocities of the 
neutrinos at Zj are frozen in time. This should improve the estimate of the matter power 
spectrum in two ways: First, by including the perturbed neutrino component the power 
spectrum will increase. Second, by switching off gravity, artificial small-scale structures 
will not form. With this method the error was decreased, especially on large scales in the 
J2m u = 0.6 eV case for Zj = 4, but the error did not approach percent level precision. 
On the contrary, in the Y2 m f = 0.15 eV case and for Zj = 49 the error increased to the 
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10-15% level. The reason is that at such a high redshift the neutrino gravitational flow 
velocities and density perturbations are very small effectively freezing the neutrinos on 
a regular grid. Because of neutrino iV-body particle shot-noise this stationary regular 
grid deflects the CDM particles and therefore reduces the amount of structure forming. 

5. Discussion and Conclusions 

We have performed a precise calculation of the effect of including neutrino dark matter 
on the matter power spectrum. The effect of thermal neutrino motion has been 
included directly, and we have shown that this effect changes the matter power spectrum 
significantly. 

Specifically, we find that the suppression of power due to the presence of massive 
neutrinos is increased by non-linear effects. Whereas in linear theory the suppression 
of power on small scales is given roughly by AP/P ~ — 8f2„/f2 m , the full non-linear 
calculation gives AP/P ~ — 9.8f2„/f2 m at a scale of k ~ 0.5 — 1 /iMpc -1 , i.e. an increase 
of about 20%. 

On smaller scales the non-linear contribution to the suppression decreases again. 
This effect has previously been noted in semi-analytic studies such as |40| . and occurs 
on a scale which depends upon the amount of CDM and the neutrino free-streaming 
length. 

The increased suppression due to non-linear effects is highly relevant for future 
high precision large-scale structure and weak lensing surveys. Even for neutrino masses 
approaching the lower bound, found from oscillation experiments, it is large enough 
to bias the estimate of other cosmological parameters. Conversely, it provides a very 
distinct signature which could allow for the detection and measurement of even very 
low-mass neutrinos. 

For k > 0.1/iMpc -1 , which include the modes for which A-body simulations are 
needed, it is a better approximation to the "true" power spectrum to neglect the 
perturbed neutrino component than to include it without thermal velocities. In both 
cases the approximation to the "true" power spectrum is better than if a standard 
ACDM model was assumed. 

But only in our lowest mass neutrino case, m v — 0.15 eV, is the error on the non- 
linear matter power spectrum made by neglecting the perturbed neutrino component 
in the iV-body simulation at the desired 1% level for the relevant modes, even for an 
iV-body simulation with a starting redshift as high as 49. The error would be decreased 
further for smaller neutrino masses. 

An alternative method for implementing the physics of neutrinos in A^-body 
simulations is to represent the neutrino component as a fluid and solve the corresponding 
fluid equations on a grid. This method would be particularly useful for neutrino masses 
close to the lower observational bound. 
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